Librerias

Primero incorporamos las librerias necesarias para poder ejecutar el presente notebook.

# Librerias
#install.packages('gganimate')
#install.packages("heatmaply")
#install.packages('gifski')
#install.packages('png')
#install.packages('units')
#install.packages('sf')
#install.packages('transformr')
require(tidyverse)
require(geosphere)
#<<<<<<< Updated upstream
require(rvest)
require(gganimate)
require(lubridate)
require(heatmaply)
require(transformr)
require(gifski)
require(png)
require(ggplot2)
#=======
require(rvest) # Cargamos el paquete
require(lubridate)
#>>>>>>> Stashed changes

Importar datos

Posteriormente importamos los datos que queremos analizar y los formateamos.

# Llamamos a los DataSets

datos_2021 = read.csv('../TP_LaboDeDatos/Data/sna_abril_2021_fixed_encoding.csv', encoding = 'UTF-8', sep = ',')

# Concateno todos los datos
datos2= read_csv2('../TP_LaboDeDatos/Data/202109-informe-ministerio.csv')

datos2$Fecha_completa = strptime(paste(datos2$Fecha, datos2$`Hora UTC`), format = "%d/%m/%Y %H:%M:%S")

En esa celda importamos los datos de los aeropuertos.

# Utilizamos la tabla que se encuentra en 'https://en.wikipedia.org/wiki/List_of_airports_in_Argentina'
# para acceder a las variables de ciudad, provincia y coordenadas de cada aeropuerto 

aeropuertos_wiki = read_html('https://en.wikipedia.org/wiki/List_of_airports_in_Argentina')
elemento_tabla   = html_element(aeropuertos_wiki,'.wikitable')
aeropuertos      = html_table(elemento_tabla)

Particularmente, necesitamos un procesamiento especial para obtener las coordenadas de cada aeropuerto en un formato que nos sea util.

# Corrijo la columna de coordenadas 

separar     = strsplit((aeropuertos$Coordinates), '/') # Divide a los strings en los lugares donde haya '/' 
coordenadas = sapply(separar, function(x) x[3])        # Me quedo solo con el 3er tipo de coordenada
coordenadas = gsub('[^0-9,.,-]','', coordenadas)           # Elimino los caracteres que no quiero utilizar

aeropuertos = aeropuertos %>% 
  mutate(lat = as.numeric(substr(coordenadas, 1, 9)), long = as.numeric(substr(coordenadas, 10, 18))) # Separo a mano latitud y longitud (revizar si esta todo en orden)

aeropuertos = filter(aeropuertos, nchar(ICAO)>1)
aeropuertos = aeropuertos[order(aeropuertos$ICAO),]

Procesar datos

En la siguiente celda lo que hacemos es a cada vuelo le agregamos los datos relativos a sus aeropuertos de origen y destino (si los tenemos). Además, calculamos con ellos la distancia recorrida por el vuelo.

for(i in 1:length(aeropuertos$ICAO)){
  inds = datos2$Aeropuerto==aeropuertos$IATA[i]
  datos2[inds,c('ciudad_origen','provincia_origen','lat_origen','long_origen')] = aeropuertos[i,c("City served","Province","lat","long")]
  
  inds = datos2$`Origen / Destino`==aeropuertos$IATA[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = aeropuertos[i,c("City served","Province","lat","long")]
}
  for(i in 1:length(datos_2021$ana)){
  inds = datos2$Aeropuerto==datos_2021$ana[i]
  datos2[inds,c('ciudad_origen','provincia_origen','lat_origen','long_origen')] = datos_2021[i,c("nam","cpr","x","y")]
  
  
  inds = datos2$`Origen / Destino`==datos_2021$ana[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = datos_2021[i,c("nam","cpr","x","y")]
  
  inds = datos2$Aeropuerto==datos_2021$iko[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = datos_2021[i,c("nam","cpr","x","y")]
  
  inds = datos2$`Origen / Destino`==datos_2021$iko[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = datos_2021[i,c("nam","cpr","x","y")]
}

datos2 = drop_na(datos2)
datos2$distancia = distHaversine(datos2[,c("long_origen","lat_origen")],datos2[,c("long_destino","lat_destino")])/1000

datos2 = datos2 %>% 
  filter(`Aerolinea Nombre` != 0)

Dividimos los eventos del dataset de vuelos en aterrizajes y despegues.

datos2$fecha_hora = strptime(paste(datos2$Fecha, datos2$`Hora UTC`), format = "%d/%m/%Y %H:%M:%S")
despegues = datos2[datos2$`Tipo de Movimiento` == 'Despegue',]
aterrizajes = datos2[datos2$`Tipo de Movimiento` == 'Aterrizaje',]

Ahora tratamos de unir los aterrizajes y despegues con el objetivo de poder identificar todo el recorrido de un vuelo. (particularmente la hora de salida y de llegada, que son los datos que no nos dan directamente en el dataset)

# dividimos en despegues y aterrizajes


matched = left_join(despegues, aterrizajes, by= c("Aeropuerto" = "Origen / Destino", "Origen / Destino" = "Aeropuerto", "Aerolinea Nombre" = "Aerolinea Nombre", "Aeronave" = "Aeronave")) %>% 
  mutate(tdif = as.numeric(fecha_hora.y - fecha_hora.x, units='hours')) %>%
  group_by(Aeropuerto, fecha_hora.x, Aeronave, `Aerolinea Nombre`) %>%
  filter(tdif > 0) %>% 
  filter(tdif < 5) %>% 
  filter(tdif == min(tdif))

matched = matched %>%
  mutate(vel_media = distancia.y/tdif)

Preguntas

Concentración del mercado

Las primeras preguntas que nos hemos planteado son relativas al grado de concentración que tienen los vuelos.

Por esto, la primera pregunta que nos interesa responder es ver cuales aeropuertos se conectan más y cuales menos. Además, notando que los casos donde el vuelo despega y aterriza en el mismo lugar son muy posiblemente vuelos de entrenamiento.

aers = unique(datos2$Aeropuerto)
mat = matrix(nrow = length(aers), ncol = length(aers))
colnames(mat) = aers
row.names(mat) = aers
for(i in 1:length(aers)){
  for(j in 1:length(aers)){
    mat[i,j]  =  sum((datos2$Aeropuerto==aers[i]) & (datos2$`Origen / Destino`==aers[j]))
  }
}
heatmaply(mat, Rowv = NA, Colv = NA,  col_dend_up=FALSE)%>% 
  layout(xaxis = list(side = "top"))

Ahora vamos a ver qué ocurre con el grafico anterior si lo pasamos a escala logaritmica (agregando 1 para evitar los -inf)

mat2 = log(mat+1)
heatmaply(mat2, Rowv = NA, Colv = NA,  col_dend_up=FALSE)%>% 
  layout(xaxis = list(side = "top"))

Como se puede ver, hay mayor cantidad de relaciones que se vuelven visibles y diferenciables gracias al cambio de escala.

Resalta el aeropuerto FDO como un importante centro de entrenamiento. Por otro lado, se observa como la conexión más fuerte es claramente entre Aeroparque y Bariloche.

La siguiente pregunta será relativa a la concentración entre distintas aerolineas.

aerolineas = unique(matched$`Aerolinea Nombre`)
apps = matrix(nrow = 1, ncol = length(aerolineas))
colnames(apps) = aerolineas
for(ar in aerolineas){
  apps[1,ar] = sum((matched$`Aerolinea Nombre`==ar))
}
apps = apps[1,order(apps, decreasing=T)]
#barplot(apps, las = 2, cex.names=0.4)
ejey <- list(title="Aerolineas", categoryorder = "array",
              categoryarray = rev(names(apps)), tickfont = list(size=4))
dfapps = data.frame(names(apps),apps)
dfapps = rename(dfapps,Aerolineas = names.apps., Cantidad = apps)
fig <- plot_ly(dfapps, x = ~Cantidad, y = ~Aerolineas, type = "bar") %>% 
 layout(yaxis = ejey)

fig

Y ahora probamos de repetir el ejercicio de expresar las distintas cantidades en escala logaritmica para ver el efecto.

appslog = log(apps+1)
#barplot(apps, las = 2, cex.names=0.4)
ejey <- list(title="Aerolineas", categoryorder = "array",
              categoryarray = rev(names(apps)), tickfont = list(size=4))
dfapps = data.frame(names(appslog),appslog)
dfapps = rename(dfapps,Aerolineas = names.appslog., Cantidad = appslog)
fig <- plot_ly(dfapps, x = ~Cantidad, y = ~Aerolineas, type = "bar") %>% 
 layout(yaxis = ejey)
fig

En este nuevo grafico se confirma la idea de que hay una enorme concentración de los vuelos en pocas empresas.

En este tercer grafico solo mostramos las más importantes para que sea entendible la información, además de que tiene una función más sintetica que los anteriores.

x = 8
apps2 = apps[1:x]
apps2[x+1] = sum(apps[x:length(apps)])
names(apps2)[x+1] = "OTRAS"
names(apps2) = paste(names(apps2),apps2)

fig <- plot_ly(data = as.data.frame(apps2), labels = names(apps2), values = apps2, type = 'pie', sort=FALSE)
fig <- fig %>% layout(title = 'Cantidad de vuelos por aerolinea',
         xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE)) 
fig

Tiempo en el aire

Las siguientes preguntas que nos hemos planteado feron las relativas a cuando los aviones estaban volando y cuando no.

La primera fue decir exactamente en qué horas del año estuvieron más y menos aviones en el aire. Para eso primero hay que agregar una hora universal a los vuelos

inicio = as.Date("2021-01-01")
matched$despegueUniversal = as.integer(-difftime(inicio,matched$Fecha_completa.x, units="hours"))
matched$llegadaUniversal = as.integer(-difftime(inicio,matched$Fecha_completa.y, units="hours")+0.9999)
fin = max(matched$llegadaUniversal)+1
histo = matrix(nrow=1,ncol=fin)
histo[,c(1:fin)] = 0
for(i in 1:length(matched$despegueUniversal)){
  inds = matched$despegueUniversal[i]:matched$llegadaUniversal[i]
  histo[1,inds] = histo[1,inds]+1
  #print(histo[1,inds])
  #for(j in matched$despegueUniversal[i]:matched$llegadaUniversal[i]){
  #  histo[1,j] = histo[1,j]+1
  #}
}
5%%2
## [1] 1
#plot(1:length(histo[1,]),histo[1,], type="s")
ts= paste((inicio + 1:fin/24), as.character(1:fin%%24))
eje <- list(title="Hora", categoryorder = "array",
              categoryarray = as.character(ts), tickfont = list(size=5))
print(typeof(eje$categoryarray))
## [1] "character"
temp = data.frame(histo[1,],Fecha =eje$categoryarray)
temp = rename(temp,Cantidad = histo.1...)
fig <- plot_ly(temp, y = ~Cantidad, x = ~Fecha, type = "bar") %>% 
 layout(xaxis = eje)

fig

Para evitar el ruido y las pequeñas variaciones diarias y semanales, vamos a repetir el grafico suavizado por día y por semana.

#plot(1:length(histo[1,]),histo[1,], type="s")

histo2 = matrix(nrow = 1, ncol = length(histo[1,])-23)
length(histo2)
## [1] 6533
for(i in 24:(fin-23)){
  #print((i-23):i)
  histo2[1,i-23] = sum(histo[1,(i-23):i])/24
}
ts= paste((inicio + 24:fin/24), as.character(24:fin%%24))
eje <- list(title="Hora", categoryorder = "array",
              categoryarray = as.character(ts), tickfont = list(size=5))
print(typeof(eje$categoryarray))
## [1] "character"
temp = data.frame(histo2[1,],Fecha =eje$categoryarray)
temp = rename(temp,Cantidad = histo2.1...)
fig <- plot_ly(temp, y = ~Cantidad, x = ~Fecha,  type = 'scatter', mode = 'lines') %>% 
 layout(xaxis = eje)

fig
#plot(1:length(histo[1,]),histo[1,], type="s")
lim = 24*7
histo3 = matrix(nrow = 1, ncol = length(histo[1,])-lim+1)
for(i in lim:(fin-lim+1)){
  #print((i-23):i)
  histo3[1,i-lim+1] = sum(histo[1,(i-lim+1):i])/lim
}
ts= paste((inicio + lim:fin/24), as.character(lim:fin%%24))
eje <- list(title="Hora", categoryorder = "array",
              categoryarray = as.character(ts), tickfont = list(size=5))
print(typeof(eje$categoryarray))
## [1] "character"
temp = data.frame(histo3[1,],Fecha =eje$categoryarray)
temp = rename(temp,Cantidad = histo3.1...)
fig <- plot_ly(temp, y = ~Cantidad, x = ~Fecha,  type = 'scatter', mode = 'lines') %>% 
 layout(xaxis = eje)

fig

En el grafico con el promedio semanal movil se puede apreciar especialmente la caida en la cantidad de aviones en el aire durante la mitad del año

Ahora cual es el horario en el que los pasajeros viajan mas

require(gganimate)
horario_vuelo = datos2 %>%
  filter(datos2$Pasajeros > 0) %>% 
  group_by(horario = hour(`Hora UTC`)) %>% 
  summarise(pasajeros = mean(Pasajeros))

table(horario_vuelo$horario)
## 
##  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
plot = horario_vuelo %>% 
  ggplot(aes(x = horario, y = pasajeros, color = pasajeros)) + geom_line() +
    geom_point() +
  transition_reveal(horario)+labs(title = 'Cantidad promedio de pasajeros a lo largo del dia', x = 'Hora', y = 'Cantidad promedio de pasajeros')+ shadow_mark()+ scale_x_continuous(breaks = round(seq(min(horario_vuelo$horario), max(horario_vuelo$horario), by = 1),1))
animate(plot, renderer = gifski_renderer(loop = T))

plot <- horario_vuelo %>% 
  ggplot(aes(x = horario, y = pasajeros, fill = pasajeros)) + geom_bar(stat='identity')+transition_states(horario,100,5)+shadow_mark()+labs(title = 'Cantidad promedio de pasajeros a lo largo del dia', x = 'Hora', y = 'Cantidad promedio de pasajeros')+ scale_x_continuous(breaks = round(seq(min(horario_vuelo$horario), max(horario_vuelo$horario), by = 1),1))
animate(plot, renderer = gifski_renderer(loop = T))

Cuales son los 3 modelos de avion mas usado por aerolinea

aerolineas_data= datos2 %>% 
  group_by(`Aerolinea Nombre`,Aeronave) %>%
  summarise(cant = n()) %>% 
  filter(Aeronave != 0) %>% 
  arrange(desc(cant))
## `summarise()` has grouped output by 'Aerolinea Nombre'. You can override using the `.groups` argument.